#include "allincludes.h"


int OneLine::LoadOverlapCoefficient(string fileName, int newBinsLengthX)
{
	//!!haven't considered about the INDEX of 15m....
	if(fileName=="")
	{
		if(newBinsLengthX==1)
		{
			fileName="OverlapCoefficient_7.5m.ovlp";
		}
		else if(newBinsLengthX==2)
		{
			fileName="OverlapCoefficient_15m.ovlp";
		}
		else
		{
			cout<<"Error : Cannot load proper Overlap Coefficient File! "<<endl;
			exit(FAILED);
		}
		fileName=WHERE_IS_INI+fileName;
	}
	cout<<"  Loading Overlap Coefficient File: "<<fileName<<" ..."<<endl;
	
	ifstream fileStream(fileName.c_str());
	if(!fileStream)
	{
		cout<<"Error : Cannot load Overlap Coefficient file!"<<endl;
		exit(FAILED);
	}
	else
	{
		fileStream.seekg(0,ios::beg);
		char buffer[1024]; //for the first 3 comment lines, do not make them too long:)
		for(int i=0;i<3;i++)
		{
			fileStream.getline(buffer, 1024);
		}
		int dataNum, fileBinsLengthX;
		fileStream>>dataNum>>fileBinsLengthX;
		//Make this clear later, about how to make the 15m resolution overlap data
		if(fileBinsLengthX!=this->binsLengthX)
		{
			cout<<"Error : BinsLength not match!"<<endl;
			exit(FAILED);
		}
		int index;
		float dist, coef;
		for(int i=0;i<dataNum;i++)
		{
			fileStream>>index>>dist>>coef;
			if(index<MAXBINSNUM && index>= 0)
				overlapCoefficient[index]=coef;
		}
		overlapEndIndex=index;
		fileStream.close();
		return SUCCEED;
	}
	return SUCCEED;  //this is really unneccessary, but the pgCC always warns about it...
}
int OneLine::LoadAfterpulse(string fileName, int newBinsLengthX)
{
		//!!haven't considered about the INDEX of 15m....
	if(fileName=="")
	{
		if(newBinsLengthX==1)
		{
			fileName="Afterpulse_7.5m.aftp";
		}
		else if(newBinsLengthX==2)
		{
			fileName="Afterpulse_15m.aftp";
		}
		else
		{
			cout<<"Error : Cannot load proper Afterpulse File! "<<endl;
			exit(FAILED);
		}
		fileName=WHERE_IS_INI+fileName;
	}
	
	ifstream fileStream(fileName.c_str());
	cout<<"  Loading Afterpulse File: "<<fileName<<" ..."<<endl;

	if(!fileStream)
	{
		cout<<"Error : Cannot load Afterpulse File!"<<endl;
		exit(FAILED);
	}
	else
	{
		fileStream.seekg(0,ios::beg);
		char buffer[1024]; //for the first 3 comment lines, do not make them too long:)
		for(int i=0;i<3;i++)
		{
			fileStream.getline(buffer, 1024);
		}
		int dataNum, fileBinsLengthX;
		fileStream>>dataNum>>fileBinsLengthX;
		//Make this clear later, about how to make the 15m resolution data
		if(fileBinsLengthX!=this->binsLengthX)
		{
			cout<<"Error : BinsLength not match!"<<endl;
			exit(FAILED);
		}
		int index;
		float dist, coef;
		for(int i=0;i<dataNum;i++)
		{
			fileStream>>index>>dist>>coef;
			if( index < MAXBINSNUM && index >= 0 )
				afterpulse[index]=coef;
		}
		afterpulseEndIndex=index;
		fileStream.close();
		return SUCCEED;
	}
	return SUCCEED;  // No branch will reach this...
}
int OneLine::LoadAtmosProfile(int ProfileNo, string fileName, int newBinsLengthX) //
{	
	//Maybe it's good to use STL vector here
	vector <DistDataSet> zTP ; //height - T - P
	stringstream zstream;
	stringstream Tstream;
	stringstream Pstream;
	if(ProfileNo==999) //Load from File
	{
		cout<<"  Loading User Defined Atmosphere Profile... "<<endl;
		if(fileName=="")
		{
			fileName=WHERE_IS_INI+"AtmosProfile.aprf";
			
		}
		ifstream fileStream(fileName.c_str());
		if(!fileStream)
		{
			cout<<"Error : Can't open atmos profile file: "<<fileName<<endl;
			exit(FAILED);
		}
		else
		{
			fileStream.seekg(0,ios::beg);
			char buffer[1024]; //for the first 3 comment lines, do not make them too long:)
			for(int i=0;i<3;i++)
			{
				fileStream.getline(buffer, 1024);
			}

			int totalNum;
			fileStream>>totalNum;

			float z, T, P;
			string tempzstring="";
			string tempTstring="";
			string tempPstring="";
			for(int count =0 ;count<=totalNum && (!fileStream.eof()); count++)
			{
				fileStream>>z>>T>>P;
				zstream<<z<<" ";
				Tstream<<T<<" ";
				Pstream<<P<<" ";
			}
			fileStream.close();
		}
	}
	else           //the following data are from LOWTRAN
	{
		zstream<<"0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 \
			 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 25.0 27.5 30.0 32.5 35.0 \
			 37.5 40.0 42.5 45.0 47.5 50.0 55.0 60.0 65.0 70.0 75.0 80.0 85.0 90.0 95.0 \
			100.0 105.0 110.0 115.0 120.0";

		switch(ProfileNo)
		{
		case 1:    //Tropical Atmosphere
			cout<<"  Loading Profile 1 : Tropical Atmosphere"<<endl;
			Pstream<<"1.013E+03 9.040E+02 8.050E+02 7.150E+02 6.330E+02 \
				5.590E+02 4.920E+02 4.320E+02 3.780E+02 3.290E+02 \
				2.860E+02 2.470E+02 2.130E+02 1.820E+02 1.560E+02 \
				1.320E+02 1.110E+02 9.370E+01 7.890E+01 6.660E+01 \
				5.650E+01 4.800E+01 4.090E+01 3.500E+01 3.000E+01 \
				2.570E+01 1.763E+01 1.220E+01 8.520E+00 6.000E+00 \
				4.260E+00 3.050E+00 2.200E+00 1.590E+00 1.160E+00 \
				8.540E-01 4.560E-01 2.390E-01 1.210E-01 5.800E-02 \
				2.600E-02 1.100E-02 4.400E-03 1.720E-03 6.880E-04 \
				2.890E-04 1.300E-04 6.470E-05 3.600E-05 2.250E-05";
			Tstream<<"299.70 293.70 287.70 283.70 277.00 \
				270.30 263.60 257.00 250.30 243.60 \
				237.00 230.10 223.60 217.00 210.30 \
				203.70 197.00 194.80 198.80 202.70 \
				206.70 210.70 214.60 217.00 219.20 \
				221.40 227.00 232.30 237.70 243.10 \
				248.50 254.00 259.40 264.80 269.60 \
				270.20 263.40 253.10 236.00 218.90 \
				201.80 184.80 177.10 177.00 184.30 \
				190.70 212.00 241.60 299.70 380.00";
			break;
		case 2:		//MidLatitude Summer
			cout<<"  Loading Profile 2 : MidLatitude Summer"<<endl;
			Pstream<<"1.013E+03 9.020E+02 8.020E+02 7.100E+02 6.280E+02 \
				5.540E+02 4.870E+02 4.260E+02 3.720E+02 3.240E+02 \
				2.810E+02 2.430E+02 2.090E+02 1.790E+02 1.530E+02 \
				1.300E+02 1.110E+02 9.500E+01 8.120E+01 6.950E+01 \
				5.950E+01 5.100E+01 4.370E+01 3.760E+01 3.220E+01 \
				2.770E+01 1.907E+01 1.320E+01 9.300E+00 6.520E+00 \
				4.640E+00 3.330E+00 2.410E+00 1.760E+00 1.290E+00 \
				9.510E-01 5.150E-01 2.720E-01 1.390E-01 6.700E-02 \
				3.000E-02 1.200E-02 4.480E-03 1.640E-03 6.250E-04 \
				2.580E-04 1.170E-04 6.110E-05 3.560E-05 2.270E-05";
			Tstream<<"294.20 289.70 285.20 279.20 273.20 \
				267.20 261.20 254.70 248.20 241.70 \
				235.30 228.80 222.30 215.80 215.70 \
				215.70 215.70 215.70 216.80 217.90 \
				219.20 220.40 221.60 222.80 223.90 \
				225.10 228.45 233.70 239.00 245.20 \
				251.30 257.50 263.70 269.90 275.20 \
				275.70 269.30 257.10 240.10 218.10 \
				196.10 174.10 165.10 165.00 178.30 \
				190.50 222.20 262.40 316.80 380.00";
			break;
		case 3:		//MidLatitude Winter
			cout<<"  Loading Profile 3 : MidLatitude Winter"<<endl;
			Pstream<<"1.018E+03 8.973E+02 7.897E+02 6.938E+02 6.081E+02 \
				5.313E+02 4.627E+02 4.016E+02 3.473E+02 2.993E+02 \
				2.568E+02 2.199E+02 1.882E+02 1.611E+02 1.378E+02 \
				1.178E+02 1.007E+02 8.610E+01 7.360E+01 6.280E+01 \
				5.370E+01 4.580E+01 3.910E+01 3.340E+01 2.860E+01 \
				2.440E+01 1.646E+01 1.110E+01 7.560E+00 5.180E+00 \
				3.600E+00 2.530E+00 1.800E+00 1.290E+00 9.400E-01 \
				6.830E-01 3.620E-01 1.880E-01 9.500E-02 4.700E-02 \
				2.220E-02 1.030E-02 4.560E-03 1.980E-03 8.770E-04 \
				4.074E-04 2.000E-04 1.057E-04 5.980E-05 3.600E-05";
			Tstream<<"272.20 268.70 265.20 261.70 255.70 \
				249.70 243.70 237.70 231.70 225.70 \
				219.70 219.20 218.70 218.20 217.70 \
				217.20 216.70 216.20 215.70 215.20 \
				215.20 215.20 215.20 215.20 215.20 \
				215.20 215.50 217.40 220.40 227.90 \
				235.50 243.20 250.80 258.50 265.10 \
				265.70 260.60 250.80 240.90 230.70 \
				220.40 210.10 199.80 199.50 208.30 \
				218.60 237.10 259.50 293.00 333.00";
			break;
		case 4:		//SubArctic Summer
			cout<<"  Loading Profile 4 : SubArctic Summer"<<endl;
			Pstream<<"1.010E+03 8.960E+02 7.929E+02 7.000E+02 6.160E+02 \
				5.410E+02 4.740E+02 4.130E+02 3.590E+02 3.108E+02 \
				2.677E+02 2.300E+02 1.977E+02 1.700E+02 1.460E+02 \
				1.260E+02 1.080E+02 9.280E+01 7.980E+01 6.860E+01 \
				5.900E+01 5.070E+01 4.360E+01 3.750E+01 3.228E+01 \
				2.780E+01 1.923E+01 1.340E+01 9.400E+00 6.610E+00 \
				4.720E+00 3.400E+00 2.480E+00 1.820E+00 1.340E+00 \
				9.870E-01 5.370E-01 2.880E-01 1.470E-01 7.100E-02 \
				3.200E-02 1.250E-02 4.510E-03 1.610E-03 6.060E-04 \
				2.480E-04 1.130E-04 6.000E-05 3.540E-05 2.260E-05";
			Tstream<<"287.20 281.70 276.30 270.90 265.50 \
				260.10 253.10 246.10 239.20 232.20 \
				225.20 225.20 225.20 225.20 225.20 \
				225.20 225.20 225.20 225.20 225.20 \
				225.20 225.20 225.20 225.20 226.60 \
				228.10 231.00 235.10 240.00 247.20 \
				254.60 262.10 269.50 273.60 276.20 \
				277.20 274.00 262.70 239.70 216.60 \
				193.60 170.60 161.70 161.60 176.80 \
				190.40 226.00 270.10 322.70 380.00";
			break;
		case 5:		//SubArctic Winter
			cout<<"  Loading Profile 5 : SubArctic Winter"<<endl;
			Pstream<<"1.013E+03 8.878E+02 7.775E+02 6.798E+02 5.932E+02 \
				5.158E+02 4.467E+02 3.853E+02 3.308E+02 2.829E+02 \
				2.418E+02 2.067E+02 1.766E+02 1.510E+02 1.291E+02 \
				1.103E+02 9.431E+01 8.058E+01 6.882E+01 5.875E+01 \
				5.014E+01 4.277E+01 3.647E+01 3.109E+01 2.649E+01 \
				2.256E+01 1.513E+01 1.020E+01 6.910E+00 4.701E+00 \
				3.230E+00 2.243E+00 1.570E+00 1.113E+00 7.900E-01 \
				5.719E-01 2.990E-01 1.550E-01 7.900E-02 4.000E-02 \
				2.000E-02 9.660E-03 4.500E-03 2.022E-03 9.070E-04 \
				4.230E-04 2.070E-04 1.080E-04 6.000E-05 3.590E-05";
			Tstream<<"257.20 259.10 255.90 252.70 247.70 \
				240.90 234.10 227.30 220.60 217.20 \
				217.20 217.20 217.20 217.20 217.20 \
				217.20 216.60 216.00 215.40 214.80 \
				214.20 213.60 213.00 212.40 211.80 \
				211.20 213.60 216.00 218.50 222.30 \
				228.50 234.70 240.80 247.00 253.20 \
				259.30 259.10 250.90 248.40 245.40 \
				234.70 223.90 213.10 202.30 211.00 \
				218.50 234.00 252.60 288.50 333.00";
			break;
		default:	//1976 U.S. Standard Atmosphere
			cout<<"  Loading Profile 6 : 1976 U.S. Standard Atmosphere"<<endl;
			Pstream<<"1.013E+03 8.988E+02 7.950E+02 7.012E+02 6.166E+02 \
				5.405E+02 4.722E+02 4.111E+02 3.565E+02 3.080E+02 \
				2.650E+02 2.270E+02 1.940E+02 1.658E+02 1.417E+02 \
				1.211E+02 1.035E+02 8.850E+01 7.565E+01 6.467E+01 \
				5.529E+01 4.729E+01 4.047E+01 3.467E+01 2.972E+01 \
				2.549E+01 1.743E+01 1.197E+01 8.010E+00 5.746E+00 \
				4.150E+00 2.871E+00 2.060E+00 1.491E+00 1.090E+00 \
				7.978E-01 4.250E-01 2.190E-01 1.090E-01 5.220E-02 \
				2.400E-02 1.050E-02 4.460E-03 1.840E-03 7.600E-04 \
				3.200E-04 1.450E-04 7.100E-05 4.010E-05 2.540E-05";
			Tstream<<"288.20 281.70 275.20 268.70 262.20 \
				255.70 249.20 242.70 236.20 229.70 \
				223.30 216.80 216.70 216.70 216.70 \
				216.70 216.70 216.70 216.70 216.70 \
				216.70 217.60 218.60 219.60 220.60 \
				221.60 224.00 226.50 230.00 236.50 \
				242.90 250.40 257.30 264.20 270.60 \
				270.70 260.80 247.00 233.30 219.60 \
				208.40 198.60 188.90 186.90 188.40 \
				195.10 208.80 240.00 300.00 360.00";
		}
	}
//	float tempDist, tempT, tempP;
	DistDataSet tempSet;
	while(!zstream.eof())
	{
		zstream>>tempSet.dist;
		Tstream>>tempSet.elem1;
		Pstream>>tempSet.elem2;
		zTP.push_back(tempSet);		//if something goes wrong here, I'd change the struct into a class with copy construct function
	}
	if(zTP.size()<2)
	{
		cout<<"Error : Can't make atmos profile with these data!"<<endl;
		exit(FAILED);
	}
	vector<DistDataSet>::iterator zTPsweeper=zTP.begin();
	float indexDist;
	float interpT, interpP;
	for(int index= FIRSTDATAINDEX;index<MAXBINSNUM;)
	{
		indexDist=IndexToDistance(index);
		if(zTPsweeper < zTP.end()-1)
		{
			if((indexDist >= zTPsweeper->dist)&&(indexDist < (zTPsweeper+1)->dist))
			{
				interpT=interp(zTPsweeper->dist,zTPsweeper->elem1, (zTPsweeper+1)->dist,(zTPsweeper+1)->elem1,indexDist);
				interpP=interp(zTPsweeper->dist,zTPsweeper->elem2, (zTPsweeper+1)->dist,(zTPsweeper+1)->elem2,indexDist);
				TProfile[index]=interpT;
				PProfile[index]=interpP;
				index++;
			}
			else if(indexDist>= (zTPsweeper+1)->dist)
			{
				zTPsweeper++;
			}
		}
		else
		{
			interpT=interp((zTPsweeper-1)->dist,(zTPsweeper-1)->elem1, zTPsweeper->dist,zTPsweeper->elem1,indexDist);
			interpP=interp((zTPsweeper-1)->dist,(zTPsweeper-1)->elem2, zTPsweeper->dist,zTPsweeper->elem2,indexDist);
			TProfile[index]=interpT;
			PProfile[index]=interpP;
			index++;
		}
	}
	return SUCCEED;
}

int OneLine::CalculateBetaM()
{
	cout<<"  Calculating BetaM..."<<endl;
	//double n0=1.000272599; //(n0-1)*10e8=6432.8+2949810/(146-lambda^-2)+25540/(41-lambda^-2);
	//Edlen 1953, I doubt about it though
	double n0=6432.8+2949810/(146-pow(LAMBDA,-2))+25540/(41-pow(LAMBDA,-2));
	//BetaM(z)=PI^2/(N0^2*lambda^4) * (n0^2-1)^2 * (6+3*rho_n)/(6-7*rho_n) * N0 * T0/P0 * P(z)/T(z)
	//where N0=2.547e25 , T0=288, P0=1013.25, rho_n=0.0284, lambda=523?
	//and PI^2/N0 * T0/P0 =1.1014e-25
	for (int i=FIRSTDATAINDEX;i<MAXBINSNUM;i++)
	{
		//Well, take care of the Units here...
		betaM[i]=1.1014e-25 * pow((n0*n0-1),2)/pow(LAMBDA/1000.0,4) * (6+3*RHO_N)/(6-7*RHO_N) * PProfile[i]/TProfile[i]*1000.0;
	}
	return SUCCEED;
}
int OneLine::IntegrateBetaM()
{
	cout<<"  Integrating BetaM..."<<endl;

	integratedBetaM[FIRSTDATAINDEX]=betaM[FIRSTDATAINDEX]*0.0075*staticBinsLengthX;
	for(int i=FIRSTDATAINDEX+1;i<MAXBINSNUM;i++)
	{
		integratedBetaM[i]=integratedBetaM[i-1]+betaM[i]*0.0075*staticBinsLengthX;
	}	
	return SUCCEED;
}
int OneLine::Initialize(LidarINI newINI,int newBinsLengthX)
{
	cout<<endl<<" --- Initializing OneLine... --- "<<endl<<endl;
	binsLengthX=newBinsLengthX;        //This is rubbish? Maybe not...
	staticBinsLengthX=newBinsLengthX;  //But this one... hmm, seems I haven't used it...

	if(binsLengthX==1)
	{
		USEABLESTARTINDEX=USEABLESTARTINDEX_7_5M;
		FIRSTDATAINDEX=FIRSTDATAINDEX_7_5M;
	}
	else if(binsLengthX==2)
	{
		USEABLESTARTINDEX=USEABLESTARTINDEX_15M;
		FIRSTDATAINDEX=FIRSTDATAINDEX_15M;
	}
	else
	{
		cout<<"Error : Can't handle this resolution!"<<endl;
		exit(FAILED);
	}
	
	memcpy(&ini, &newINI, sizeof(LidarINI));
	SINALPHA=sin(ini.alpha*PI/180.0);
	COSALPHA=cos(ini.alpha*PI/180.0);

	LoadAtmosProfile(ini.atmosModel);
	CalculateBetaM();
	IntegrateBetaM();
	LoadOverlapCoefficient("",newBinsLengthX);
	LoadAfterpulse("",newBinsLengthX);
	
	return SUCCEED;
}

